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Abstract An inverse problem of wave propagation into a weakly laterally inhomogeneous 
medium occupying a half-space is considered in the acoustic approximation. The half- 
space consists of an upper layer and a semi-infinite bottom separated with an interface. 
An assumption of a weak lateral inhomogeneity means that the velocity of wave prop- 
agation and the shape of the interface depend weakly on the horizontal coordinates, 
x = (xi,X2), in comparison with the strong dependence on the vertical coordinate, z, 
giving rise to a small parameter £ << 1. Expanding the velocity in power series with 
respect to £, we obtain a recurrent system of ID inverse problems. We provide algorithms 
to solve these problems for the zero and first-order approximations. In the zero-order ap- 
proximation, the corresponding ID inverse problem is reduced to a system of non-linear 
Volterra-type integral equations. In the first-order approximation, the corresponding ID 
inverse problem is reduced to a system of coupled linear Volterra integral equations. 
These equations are used for the numerical reconstruction of the velocity in both layers 
and the interface up to 0(e 2 ). 

Key words: Acoustic equation, inverse problem, velocity reconstruction. 

1. Introduction 

The inverse problems of geo-exploration implies investigation of domains inside the earth's 
crust which contain gas, oil or other minerals as well as determination of their physical 
properties such as density, porosity, pressure, and so on. The most important problem is 
determination of the boundaries of the domains which gives information about its loca- 
tion, volume and makes possible to predict costs and outcomes in exploitation. In practice, 
in e.g. oil exploration and seismology this inverse problem is mainly solved by means of 
the map migration method. The map migration method assumes the knowledge of the 
velocity profile above the interface and provides robust, effective numerical algorithms 
which are quite stable with respect to variations of the velocity and complicated geometry 
of the interface. Therefore, it has attracted much attention of mathematicians and geo- 
physicists, see e.g. Q, lE3L lETl . 03, with new approaches continuing to appear e.g. 

However, an absence of an independent way to recover the velocity profile above the 
interface may hinder the map migration techniques. This makes crucial for geosciences 
to develop algorithms to solve an inverse problem of the velocity reconstruction from the 

measurements of the wave field on the ground surface, z = 0. 

l 
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This gives rize to a fully non-linear dynamic multidimensional inverse problem. To our 
knowledge, the only method valid for arbitrary inhomogeneous velocity profiles is the 
Boundary Control method (BC-method), see e.g. Q, EOll . However, at the moment, this 
method is developed only for smooth velocity profiles, in the absence of discontinuities. 
Moreover, existing variants of the BC-method do not possess good stability properties 
and have a rather poor numerical implementation, e.g. Q. 

On the other hand, in the case of a pure layered structure, with the properties depending 
only on depth, z, inverse problems of geophysics are often reduced to one-dimensional 
inverse problems. The theory of one-dimensional inverse problems which goes back to 
the classical results obtained by Borg, Levinson, Gel'fand-Levitan, Marchenko and Krein 
of late 40-th - 50-th is still an area of active research with new theoretical results and 
numerical algorithms appearing regularly in mathematical and geophysical literature, see 
e.g. ED. E2, ED, E3, G3, ED- 03 to mention just a few. 

Returning to geophysical applications, a pure layered structure occurs rare. However, in 
many cases in seismology and oil exploration the parameters of the earth are not described 
by arbitrary functions of 3 dimensions having jumps across arbitrary 2 dimensional sur- 
faces. Rather the properties of the medium depend mainly on the depth, z, with only slow 
dependence on horizontal coordinates, x = (x\,X2), i.e. depend on ex, £ << 1 rather than 
x (in seismology, often e ~ 0. 1 — 0.2) and we deal with a weakly laterally inhomogeneous 
medium (WLIM). The importance of WLIM is now well-understood in theoretical and 
mathematical geophysics. There are currently numerous results on the direct problem of 
the wave propagation in WLIM, see e.g. O , EH, ED, (TT), 0, E3- However, to our 
knowledge, except lfT3l. there is currently no special inversion techniques for WLIM. In 
this paper we develop and numerically implement an inversion method for WLIM based 
on different ideas than those in [13]. This method provides, for inverse problems occur- 
ring in important practical applications, a technique which inherits some practically useful 
properties of the one-dimensional inverse problems, e.g. robustness and fast convergence 
rate of numerical algorithms. 

In the acoustical approximation we use in this paper, the wave propagation in WLIM is 
described by the wave equation, 

(1) u tt -c 2 (z,x)A z ^u = 0, z>0 
with 

(2) c(z,x) = c (z) +£ < x,ci(z) > +■■■■ 

Due to ©, the response of the media, measured at z = 0, may also be decomposed in 
power series with respect to £. Analyzing this decomposition, we obtain a recurrent 
system of inverse problems for subsequent terms c p (z). The inverse problem for cq(z) 
reduces to a family of one-dimensional inverse problems. We employ as a basis to solve 
this problem the method coupled non-linear Volterra integral equations developed by 
Blagovestchenskii, e.g. |6|, which goes back to H. Note that H, |6| provide a local 
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reconstruction of co near the ground surface z = 0. In this paper, we pursue the method 
of coupled non-linear Volterra integral equations further making it global, i.e. providing 
reconstruction of co(y) from inverse data measured for t E (0,2T) up to the depth T in 
the travel-time coordinates, y (see formula (fl6l) and obtain conditional Lipschitz stabil- 
ity estimates for this procedure. The results of the zero-order reconstruction serve as a 
starting point to develop a recurrent procedure to determine, from the inverse data, further 
terms in expansion 0. Such procedure, for the inverse problem of the reconstruction 
of an unknown potential q(z,£x) rather than c(z, £*),was suggested by Blagovestchenskii 
1 5 1, |6|. The method of |6| is based on the use of polynomial moments with respect to x. 
We suggest another approach based on the Fourier transform, x = (xi,X2) — > £, = (£i > ^2), 
of the equation (HJ). This makes possible to utilize the dependence of the resulting prob- 
lems on £, in order to solve the recurrent system of ID inverse problems. The coefficients 
ci(z),C2(z), ■■■ may then be obtained as solutions to some linear Volterra integral equa- 
tions. Having said so, we note that the integral equations for higher order unknowns 
contain higher and higher order derivatives of the previously found terms, thus increasing 
ill-posedness of the inverse problem. This is hardly surprising taking into account well- 
known strong ill-posedness of multi-dimensional inverse problems JT9), ll26ll . What is, 
however, interesting is that within the model considered there is just a gradual increase of 
instability adding two derivations at each stage of the reconstruction algorithm. 

In this paper we confine ourself to the reconstruction only of co, c\. Reconstruction of 
the higher order terms c p , p > 1 is, in principle, possible using the same ideas as for c\, 
although is more technically involved. However, in practical applications in geophysics 
the measured data make possible to find inverse data only for co, c\ . Indeed, (HJ), © imply 
that, for the type of the boundary sources we consider, 

u(z,x,t) = uo(z,x,t) +eui(z,x,t) + e U2(z,x,t) + ... , 

so that the measured data at z = 0, 

u z \ z =o(x,t) =R(x,t) = R (x,t) +eRi(x,t) + e 2 R 2 (x,t) + ... . 

What is more, R p (x, t) are even with respect to x for even p and R p (x, t) are odd for odd p. 
Clearly, in real measurements R is not given as a power series with respect to £. However, 
using the fact that Rq is even, and Ri is odd, we can find Rq and R\ up to 0(e 2 ) from the 
measured R. 

The plan of the paper is as follows. In the next section we give a rigorous formulation of 
the problem and provide a general outline of the perturbation scheme for this problem in 
WLIM. In section 3 a modified method of coupled non-linear Volterra integral equations 
to reconstruct Co (z) is described. In section 4 we derive a coupled system of linear Volterra 
integral equations for ci(z). Section 5 is devoted to the global solvability of the non- 
linear system for co(z) and conditional stability estimates. In section 6 we test the method 
numerically for the two dimensional case (inhomogeneous half-plane). As we stop with 
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Co, c\, in practical applications this would result in an error of the magnitude 0(e 2 ). The 
final section is devoted to some concluding remarks. 

2. Formulation of the problem 
Let us consider the wave propagation into the inhomogeneous acoustic half-space 



(3) 



1 



n (z,£x)u tt div(pVw) = 0, x=(xi,X2), 



where n is the refractive index, n = c 1 . Above and below the interface z = h(ex) the 
refractive index is described by 

nW(z,ex), 0<z<h(ex), 



n(z,ex) = < (2) ) ; ., \ 

{ nW(z,ex), z>h(ex) 



where £ is a small parameter (0 < £ 1) which characterizes the ratio of the horizontal 
and vertical gradients of n. We assume that the density is piece- wise constant, 

f pi, 0<z<h(ex), 
\ p 2 , z>h(ex), 

with known pi, p2. As we deal with WLIM, the following notations are used, 

(4) n 2 (z,£x) =nl(z) +£ <x,n(z) > +0(£ 2 ), n(z) = (ni(z),n 2 (z)), 



(5) h(ex) = h + e < x,h > +<9(£ 2 ), h = (hi,h 2 ), 

where <, > means the scalar product. 

We assume that u = for t < 0, and the boundary condition is 



(6) 



»(*)/(*) 



z=0 



Pi 



«o(0) 



/(f) = 8(f) or /(f) =0(0, 



where 8(x) and 0(f) are 8-function and Heaviside function, correspondingly. On the in- 
terface between these two layers we assume the usual continuity conditions, 



(7) 



[u] 



du 

dn 1 



= 0, 



z=h(ex) 



z=h(ex) 

where [. . . ] stand for a jump across the interface. Using the Fourier transform with respect 
to x 



(8) 



R 2 



U can be expanded into asymptotic series, 

CO 

(9) U(z,^t)^^E n i n U^(z,^t), 



n=0 
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where all functions CA n '(z,£,f) are real. Using decompositions ©, ©, it is easily seen 
that t/( w ) (z, ^, ?) are even functions with respect to £, for even n and odd for odd n. 

Our goal is to reconstruct the refractive index and the shape of interface, namely, the 
functions no(z), «l,2(z) and constants ho, h\_2, from the response data collected during 
time < t < IT, 

du 



dz 



■R{x,t,E), 



z=0 



with 



R(x,t,e) « 



77=0 



Decomposing the wave equation (HJ) and interface conditions © with respect to £, we 
obtain initial-boundary value problems for and U^ l \ The zero-order problem is 



(10) ^(z)^ 0) -e4 0) + I^ (0) =o, |$| 2 = g+& f/ (0) 

with the interface continuity conditions 

(ID [t/ (0) ] 



z=0 



Pi 



«o(0) 



Z=/l() 



o, [V>] 



0. 



z=^o 



and inverse data of the form 



(12) 



(0) 



■ ro(t,Z,) = J cos(<^,x>)Ro(x,t)dx. 

R 2 



The first-order problem is 

nUz)U^-U^ + \^U^=<n^U t ?>, U< 



(13) 



with the interface continuity conditions 
(14) 



0. 



z=0 



' Z=h 

and inverse data of the form 



= 0, 



Uu 



(1) 



<h,V % U^>+U^<U> 







z=^o 



(15) 



(i) 



z=0 



ri(?,^) = y sin(< ^,x >)R\(x,t)dx. 

R 2 



In this paper we confine ourselves to the reconstruction of only hq, n and ho, h. The 
reconstruction of the higher order terms is, in principal, possible using the same ideas as 
for h and h, although is more technically involved. Moreover, in practical applications 
in geophysics the measured data make possible to find the inverse data only for no, h 
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and ho, h. Using the experimental data of the response R(x,t), one cannot expand it into 
power series with respect to £. However, as Rq is even and Ri is odd with respect to x, we 
have 

ro(t,Z,) = Jcos(< ^,x >)i?(jc,?)^+0(e 2 ), r\(t,Q=E~ l j sin(< ^,x >)i?(jt,f)<£x;+<9(£ 2 ). 

R 2 R 2 

So, the only quantities we may observe from the experiment are ro(^,0 up to 

error £ 2 . Thus, although the inverse problems for the higher-order terms in © can be 
considered analytically, inverse data for them are not available from the measurements. 
In this connection, we do not discuss higher approximations in this paper. Also the exact 
value for the parameter £ is not known from the experiment. However, this quantity may 
be evaluated numerically using inverse data of response R(x,t,e). For example, one of the 
ways to calculate £ is as follows 



£ = sup 



J sin(< l^x >)R(x,t)dx 

R 2 



0<t<2T, |^| <Uw, 



where %max is a positive real number of order 1. It is clear that £ is qualitatively related to 
an extent of the medium being a weakly laterally inhomogeneous one. 



3. Inverse problem in the zero-order approximation 



3.1. Half-space. In this section we describe an algorithm to solve the inverse problem in 
the zero-order approximation for an inhomogeneous half-space. Namely, we will describe 
an algorithm to determine n$(z) from the response ro(t, ^) . This is a generalization of the 
approach developed first by Blagovestchenskii @|, see also ©. The crucial point of the 
approach is the derivation of a non-linear Volterra-type system of integral equations to 
solve the inverse problem. 

Let us introduce two new independent variables 



(16) 



y 



n (z)dz, a(y) 



o 



rco(z(y)) 
P 



The function o(y) is called acoustic admittance while y is the vertical travel-time. Then, 

Let f(t) = 8(f). Then, the boundary condition and response are 

U(0) 



8(f) 



y=0 



Pi 



*o(0) 



(0) 



y=Q 



n (0) 
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Let us change the dependent variable 

(18) V i(y,0 = VW)U {0 \ V2 Cy,0 = ^ + ^, 

thus reducing the second order PDE for to a system of two PDE of the first order 
(m f ¥l + ¥iv = ¥2 

where 

(>/555)" l^l 2 



(20) q(y,%) 



2 • 



Expressions in the left-hand side of (fl9l) are the total derivatives of and \|/2 along 
corresponding characteristics (see Fig.l). Integrating the first equation along the lower 
characteristics, and the second equation - along the upper one, we obtain for t > y 



(21) 



where 



y 

Vi(y»0 = JV2{T],t + T\-y)dT], 



y 

V2(y,0 = -/4(r|,^)¥i(TV-r|+y)rfri+s(f-f)^) 




(22) ,(a)=5'(o + &, )+ ro( ^ } 



2«o(0) v/pm (0) 

is bounded as ? — > due to the cancellation of S'(j)— and 8(f)— singularities in the right- 
hand side of d2"2l. 

The system d2T1l is not complete to solve the inverse problem as it has three unknown 
functions \)/2(;y,0 and and just two equations. We use the progressive 

wave expansion, see e.g. [6|, Qui . 

t 

(23) U^(y,t) = £ /„(?-y)t/i 0) (y), / n+ i(0 = / /„(0<&, / (?) = 8(0 

to derive the third equation. The amplitudes uh°^ are found from the transport equations 
with e.g. 

(24) U® = ^L= = —}==[ q(T\)d % 
so that 

V2(y,t) = -Q(t-y)q(y) + .... 

Here and later we often skip the dependence of q and other functions on £, when our 
considerations are independent of its value. 
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Substituting the second equation of (I2TT) into this relation, and taking t — ► y + 0, we obtain 
the third desired equation. We summarize the above results in the follwing 

Proposition 3.1. Let 1(^2 (y, t) is obtained from the solution (y, t) of the wave equation 
J77t by formula rf7§l) . Then \\fi^(yj) together with q(y) satisfy the following system of 
non-linear Volterra-type equations, 



(25) 



Vi(y,0 = fV2{%t + T\-y)d% 
o 

y 

V2(y,0 = -J^(ri)x)/i(ri,f-ri+y)Jri+^(?+y), 
o 

y 

qiy) = -2jq(r\)\]fi(r\,2y-r\)dr\ + 2g(2y). 
o 



This non-linear Volterra-type system of integral equations, in the sense of Goh'berg-Krein 
HHI, allows us to determine, at least locally, the function q(y) from the response g(t) for 
< t < 2T . In section 5 we return to the question of the global solvability of system (E51 . 
Having solved this system for two different values £i and ^ we determine the refractive 
index in the zero-approximation, 



(26) 



My) 
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which must be used together with 



z 



dy 



o 



My) 

The refractive index no(y) may also be determined by solving system d25t just for a single 
value of ^. It can be done by integration the differential equation for w(y) = v/o(j) which 
follows from d20l) 

pjw 3 

The general solution of this equation may be represented in the form, see Q. 



w (y) = J L A u w i(y) w j(y)i 

V U=l,2 

where wi^Cy) make a fundamental system of the corresponding homogeneous equation. 
Constant real symmetric matrix Ay satisfies 

det||A ; - J ||W 2 (wi,w 2 ) = -^, 

Pi 

where 

W("Wi,W2) = W\w' 2 — M>\W2- 

Given w(0) and w/(0) we are then able to determine all entries of the matrix A. 

3.2. Two layers problem. Consider a wave propagation through the interface, y = y(ho) — 
L, using the singularity analysis of the incident, reflected and transmitted waves. Then, 
in the layer < y < L, the singularities of L^ ) given by the incident and reflected waves, 
are 

^(o) = 8('-y) , dyj ^ t+ y- 2L ) + 

For the transmitted wave, the WKB asymptotics is given by 

8 (,-y) | 9 ( f -y) 



where 



Qo(y) = J q(T\)d% Qiiy) = J q(i\)dr\. 
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Unknown quatities do, d\ and so, s\ are the reflection and transmission coefficients. Using 
the interface continuity conditions (fTTT) . leads to a linear system of equations for do, d\ 
and so,s\ so that 



(27) do 



so 



1 + 



^1 



a'_ ( 1 + Jo) + <2o (L) (a_ - o+) - sooV a /| 




a+ + a 

oLCl + db) + 2fio(L)o_ - jooVVS 



a_ 1 + 




o+ + a 



where a± and a' ± are the limiting values of o(y) and o'(j) at y — > L — and y — > L + 0, 
correspondingly. 




Figure 2. Characteristic lines in the case of two layers. 



Analyzing the incoming singularities at y = 0, we can find L and, therefore, h — and also 
do and d\. This will give us o+ and & + . The next step is reconstruction of the velocity 
at large depth L <y < T, which is illustrated on Fig. 2. This reconstruction requires the 
data observed for < t < 2T. At this stage we again apply the non-linear Volterra system 
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of integral equations <l25t in the new frame y',T[',t' (y = y' +L, r| = T|' + L, t = t' + L) 



(28) 



= J M[ 2 (Tl',* 7 + if - / ) di\' + (0, t'-y'), 



y 

= - / qWvi Cn', * ' -r['+ y w + y+ (o, t' + y ) 



y 

9 (y) = -2/ ? (Ti / )vi(ii , ,2y-Ti / yTi / +2 ¥ +(o,2y). 





These require knowledge of (0, — 3/) and \|/ 2 (0,f' + y) at the vertical segment AB 
as the limits \j/j~ 2 (0,?) = lirn\|/i 5 2(y,*) for y — > L + 0. To this end, we first determine 
\|/] _ 2 (0, t) = limv^i 5 2(y, t) for y — > L — using system (l2lTl . Employing the interface conti- 
nuity conditions 



¥1 _ ¥i 



+ 



,+ 



we obtain the required data, (0, t' - y') and \|/ 2 (0, t' + y') at A.B. 



4. Inverse problem in the first-order approximation 

4.1. Half-space. In this section we describe an algorithm to determine n(z) in an inho- 
mogeneous half-space from the response n(?,^). For convenience, we integrate all the 
functions of the zero-order problem with respect to time t so that 



U(0) 



y=0 v »o(0) 



Using the variables y and o(y), we rewrite the problem (fT3t . (fl4l in the form 

(29) 4" -J-M <,(y)uP) + W V m = \<n, V<> >, 

o(y)3yV / "0 n o 



U 



(i) 



y=0 



0, t/ y 



(1) 



y=0 



no(0) 



h{t£) = J n{t£)dt. 



We start with the progressive wave expansion for As by (l23l) . (l24h . 



8(?— y) /" dr\ 



<«,v 4 nf»>=-<«,5>»2=4/ 

\/a(VI >/ 



+ ... . 



the progressive wave expansion for t/*- 1 -* (y, f) has the form 

uW(y,t) = Q(t-y)A (y) + (t -y)+Ai(y) + .... 
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Therefore, equation (l29b implies, in particular, that 



(30) 



Let G(y, T|,?;£) be the space-causal Green's function, 



(31) G„ - (a(y)G^ + i^G = 8(f )8(y - r|), G = Oif y < T|. 



Then 

(3ZP>(y,y + 0;S) = 
? V <«vA T'E)> 1 r 

/rfn / dxG(y^y-x;^) ' % " \ h ' w — dxG{y,0,y-v£)h{x&) 

J J ni(r\) no(0) J 

ri uv o 

where we now write down explicitly the dependence on ^. 
Introducing new unknown functions, 

(33) (p( y ) = (91 (y) ,92 60) = n(y)p(y) 

and employing the progressive wave expansion for we see that 



<<p(y)^>=—l^c(yj J dT] J d%G(y,T\,y-X&) 



d y[ J J nm)p(T\) 



(34) 

Finally, the desired system of linear Volterra integral equations may be obtained from 
(EH) by differentiation and setting £, equal to e.g. £i = (« 5 0) and ^2 = (0, a), where a =^ 0. 
Namely, we obtain 

Proposition 4.1. Le? <p(y) = (9i(y),92(y)) £e given by Oil where the scalar factor p(y) 
depends only on the already found no(y). Then (p(y) satisfies the system of linear Volterra 
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equations 



<9(n),V^S° ) (ri,2y-Tl,y > 



a%(y) = 2 / dT[Gi(y,T[,r[-yXi) 



*o(0) 



no(t\)p(vi) 



2y-r\ 

Gi(y,0,-y,Z,dh(2y,Z,i)+ I dr\ [ dzG 2 {y,r\,y-z,\ 



2y 



<<j>(Tl),V^° ) (Tl 1 T 1 ^)> 

*oCn)/>0i) 



no(0) 



j-yjxG 2 (y,0,y-x,yn(x,y, »=1,2. 



where 



Gi(y,T\,t-x;^) = y^a(y)G(y,T\,t-x;^) G 2 (y,T\,t -%;£,) = G lt (y,T\,t -%;£,) + G ly (y,i\,t -%■£). 

Observe that, due to the first equation in (l25t . 

limV^f (ti,x)=0, 
r|^0 

so that the system of integral equtions in Proposition ^. H is not singular. 

4.2. Two layers problem. Let us consider the inverse problem for the first-order approx- 
imation in the case of a layer and a semi-infinite bottom separated by an interface. Our 
goal is to derive a coupled system of linear Volterra integral equations similar to that in 
Proposition 14. II to reconstruct n(z) and also to determine constants h characterizing this 
interface. Recall the formulation of the corresponding problem using y variable, see (IT~3l> . 
<d, 

(35) u;, " - — — i o( v ) u; " ) + ^ u ' ! ■' - — < n . v , u;; : > . 



<3{y)dy\ J Tig n 2 Q s 



[/(I) 



0, £A 



(i; 



y=0 



y=0 



no(0) ' 



U W -n <h,V^u} 0) > 



0, 



aMt/i 1 )- (pa 2 (y)<7,,V^ 0) > +pa(y)a'(y)<^V^ 0) > -tf<°)<M>) 

As before, we start with the singularity analysis of the transmitted and reflected waves. 
The progressive wave expansion of the inhomogeneous term in (133b is given by 



0. 



Z,i8(t-y) f dr[ 8(t+y-2L) (\ d 



d L M - - ■■ , : | . 



«oCn) 



+ ... 



0<y<L, L<* <2L, i = 1,2, 
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a f/ (Q) g(*-y) fi g yi E 



dr\ 



+ . 



y>L,L<t<3L. 



Then the sum of the incident and reflected waves is given by 



(36) t/ (1) 



where 



n Q(t-y) 



< nX > p(r\)dr\ + 



Q(t+y-2L) 



V a (y) 

< y < L, L < t < 2L, 



< n,pi > dr\+D +... 



i na, . 



2«§(y)V2a^ 
The transmitted wave is then 



dr\ 



0<y<L. 



(37) 
where 



< n,;?2 > dri + S +..., y>L,L<t<3L, 



2»5W 

and coefficients di , si are defined in (|2"71) . 

Substituting these expressions for into the interface continuity relations in d35t . and 
solving the corresponding system of linear equations, we find that 



(38) 



1 - ft / < n & > p(y)dy + ^ < h,Ti > + < h,T 2 > 



1 + - 



2 J < n£>p(y)dy+ < h,r 2 -Ti > 



Here the vectors Ti and T 2 are given by 

L 



rf = Pio. 



plO 



*oCn) 

*oCn) 



/or i a 

V a + 2 



c P2a + -^ l5 , = 1,2. 
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From the jump of at y = at the time t = 2L, we can find D. Now equation (l3*%b with 

4i = («,0) and ^2 = (0,a) make possible to find h, 

(39) 

L L i 



2ap 







Clearly, at this stage we can also find S. 

Let us describe an algorithm to determine n(z) in the second layer. The boundary condi- 
tions \ y =L+o and the response Uy 1 ^ \ y =L+o for the transmitted wave can be found using 
r\ (t, ^) and already known, in the upper layer, no(y) , n(y) together with the interface con- 
tinuity conditions in d35t . Using the coordinates y' ,r\'j', y = y+L,T| =T\' +L, t =t' + L, 
see Fig. 1, we have the progressive wave expansion for 

y 

= ^=f [<n&>p2{i\)di\+s) + ..., P2(y')=^P2(y'). 

Employing the Green function (l3lT) in the second layer, we obtain 

p 2 {y')<n{y')X> = 

d S rrrJ jr . , V' , , , , < ^(ii^v^^^) > 
^1^550(7*1 / ^(v,,,v-x) ^ 

^ ri' u 

2/ / a' \ V 

-1 JxG(/,0,/-x) U 2 (x) + ^±xi(x)J + / dvG^y flj -x)n{x) 



where 



+ [Giy' , 0, y )xi (0) - G(/ , 0, -y ) X i (2y )] 



Xi(t) = U% =L , x 2 (t) = u} l) \y= L 



Finally, the desired coupled system of linear Volterra integral equations determining n(y) 
below the interface may be obtained from this equation by setting e.g. = (a,0) and 
^2 = (0,a) (compare with Proposition l4.ll) . It is worth noting that these integral equations 
are not singular. 



5. Global reconstruction and stability 



As mentioned in section 3, the non-linear Volterra-type system of integral equations d25l) 
may, in principle, be solvable only locally. In this section we show that, assuming a priori 
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bounds for q(t) in C(0,T), it can be found on the whole interval (0, T) using a layer- 
stripping method based on d25t . What is more, we will prove Lipschitz stability of this 
procedure with respect to the variation of the response data g(t) in C(0, IT) . 

We first note that if | \q\ \c(o,t) < M then solving the direct problem gives 

(40) ||¥i,2||c(A r )> \k\\c(0,T)<A=A{T,M), 

where we can assume A > max( 1 , M) . Here Aj and generally A^, < b < T, is the triangle 
in R 2 bounded by the characteristics r\ = x, r\ + x = 2b and the axis r\ = 0. Assume that 
we have already found \\f\_2 in At and q in (0,b). 

Lemma 5.1. Let % > satisfies 

(41) l< ' 



16(1 +A) 2 ' 

Then system rt25l ) uniquely determines \|/i,2(y,f), q(y) for (y,t) G A b+ x, y G (0,b + X), re- 
spectively. Moreover, these functions can be found from system X25\) by Picard iterations. 

Proof Utilizing new coordinates y r =y — b,T\ , = T\—b,t , = t — b (compare with the end 
of section 3), we obtain a system of non-linear Volterra-type equations similar to (I2"%1) . 



(42) 



¥i(y,o 



JT ^2 (ti' 5 S + rf- y')dT}' + gl (/, , 
o 

y 

- / qfa'ftW OV/ - Tl 7 +/)^ 7 + #2 (/ , *') 


y 

-2/^(ri / )v)/i(ri / ,2y-ri / )rfri / + g3(y). 



Hereby/) = ¥l (0/-/), g 2 (y',t') = \|/ 2 (0,f' =j'), g 3 (yO =2^(0,2/) maybe 
found from the "direct" system (|2TT) and satisfy \gi\ < A. 

Using the first equation in (@2K we eliminate to get 
(43) 

¥ 2 (y,0 = fgi(ri',^-Ti'+y)+/ T1 V 2 (^,f'-Ti'+^ 7 l^'+g 2 (y ) o, 







-IJqtf) gi(Tl',2y-Tl')+/o T1 ¥2(^y-Tl / + ^ 7 *l' + *3C>0 
L J 



We rewrite this system in the form, 

(V2>9)=*5i(Y2,9), 

where ^ is a non-linear operator in C(A^ x C(0,A.) determined by the right-hand side 
d4*3l . Let us show that is a contraction in 

V 2 a=B 2 a(C(Ax)) x*2a(C((U)), 
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where B p (-) is the ball of radius p in the corresponding function space and we use the 
norm 

||(\|/2,4)||=max(||\|/ 2 ||, |M|). 

Considering 

K x (y 2 ,q) -K K (y 2 ,q) = (K X (y 2 ,q) -K x (y 2 ,q)) + (K x (y 2 ,q) - 
it follows from (Hoi that in D 2 a 

||^(\|f2^)-^%»^)ll<2AA,(lk-^ll+2A.|IV2-%||+2A,|k-$||). 
As X satisfies (HTb. this implies that 

\\K x (y 2 ,q) -K X (ty 2 ,q)\\ < |ll(V2, q) - 



Similar arguments show that K x : D2A — > D2A- Thus ( Hot has a unique solution in D24 
and, since d28t is a Volterra-type system, in C(A X ) x C(0, X). 

QED 

Next we prove the Lipschitz stability of the non-linear Volterra system (l25t . Let us denote 
by % 2(jj 0' ^(^ an ^ g (0 small variations of the corresponding functions. 



Lemma 5.2. Let ^1,2, q satisfy a40v . There exist cq = cq(T,A), £0 = £ — 0(T,A) such that 
for \\g\ |c(o,2r) < £ < £0 the corresponding system ( 1251) with g + g instead ofg has a unique 
solution 11/1,2 + V1.2, q + q, and 

(44) I |?l,2||c(A r ): 1 |«| lc(0,2T) < c o£- 

Proof Substituting expressions for yi :2 (y,t) +fyi, 2 (y,t), q(y) + q(y) and g(t) + g(t) into 
d25t , we obtain the corresponding system in variations 

( y 
%M = fV2(%t + T\-y)d% 


y 

¥2(j,0 = -Jq(T\)%(f\,t-T\+y)dr\ 



y 

-fq(^)(w(%t-^+y)+%('t\,t-T\+y))dT\+g(t+y), 


y y 
q(y) = -2/^(ri)\)/i(ri,2v-ri)Jri-2/(^(ri)+^(ri))v)/i(ri,2v-ri)Jri + 2|(2v). 


Let I \g\ I < £. Then, we have 

y 

\fi(y,t)\ < f\%(%t+r\-y)\dT\, 


y y 
!%()>, *)| < 2A/|tpri(Ti,J-TH-y)|rfn+A/|9(Tl)|rfn + e, 


|§(y)| < 4A/|^(ri)|Jri + 2A/|%(ri,2y-ri)|Jri + 2£. 
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Substituting the upper inequality into the second one to replace |\j/i(r|,J — T| +y)\, we 
obtain 

y 

\y 2 (y,t)\<2A y |\j/ 2 (ri,x)|JriJx + 2A j \q{y\)\dx\ +e, 

Ay,t 

where A y j is the the triangle in R 2 bounded by the characteristics % = r\ + t—y,T = t+y—r\ 
and the axis T| = 0. Similarly, 

y 

\q(y)\<4A J \q(r\)\dT\ + 4A J |\j/ 2 (r|,T)|^r|^T + 2e. 

Ay, 



Let 

Then, we have 



as 



p(f\) =max T |\j/ 2 (ri,x)|, (%x)eA T . 

y y 
\y 2 (y,t)\<4AT J p(r\)dr\ + 2Aj \q(t\)\dT\ + e, 

o o 

y y 

|«0)| <4A J \q(r\,t,)\dr] + 8AT J p(r|)Jr| + 2£, 



y y 
J p(r\)dr]dx = 2 J p(j\)(y-T\)<m<2T J p(r\)dr]. 

AyJ 

Moreover, it holds that 

y y 
p(y)<AAT J p{y\)dx\+2A J \q{r\)\dr\ +e. 
o o 

Introduce 

p(y)=p(y) + \q(y)\- 

Adding the last two inequalities for p{y) and \q{y) | yields 

y 

p(y)< J(12ATp(r))+6A\q(r))\)dr) + 3E, 
o 

or, with C = 6Amax( 1,27), 

y y n 

p(y)<C J p(ri)Jri + 3e<C j IcJ p(r]i)dr\i + 3e)dr] + 3e 



o o 



=j / p(Tl)(y-Tl)</Ti + (Cy+l)3e. 
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Continuing this process, we come to the estimate 

y 



C n f ( C" 

P(y)<— - / p(r[)(y-T\) n dT\+iy n —+y 



»-._ + ... + i l 3 ! . 



Let now po = maxp(y) for y G [0, T] . Then, 

■<n-\ 



s-<n+ 1 y n+ 1 

p(y)<Po , * +3 g %, 



i.e. 

(CT) n+l „ rr 
Po<P0 7 ' , +3e cr e. 
(n+1)! 

Clearly, for sufficiently large n it holds that 

(CT) n+1 1 
(n + 1)! " 2' 

and then, 

Po < 6<? C7 £. 
Thus, assuming | \g\ \ < £, we obtain that 

(45) 1 1^| | < 6e cr £, ||\j/ 2 || < 6e cr £, 1 1% 1 1 < 67e cr £, 

which implies (14*41) with 

c = 6max(l,7> CT , £ = ^— — e" CT (A -max(||\|/ 1;2 ||, \\q\\)). 

o max i , i j 

QED 

Observe that when g = 0(1) estimate (l45t does not imply, for sufficiently large T, that 
\j>i,2 + Vi,2, <? + <? satisfy (Hull . This is a manifestation of a possible non-convergence of 
the Picard method for large g. 

Observe also that the conditional stability, in C(0, T), of the inversion method based on 
the Volterra-type system (1251) implies the conditional stability, in C 2 (0, T), of the original 
inverse problem of the reconstruction of no(y). Indeed, if no(y) is a priori bounded in 
C 2 (0, T), i.e. q(y, ^) is a priori bounded in C 2 (0, T) for bounded ^, the inverse map, 
g(t, ^) — >• no(y), is Lipschitz stable from C(0, 2T) to C(0, T). Due to (|20T) . this implies the 
Lipschitz stability of the above map from C(0, 27") to C 2 (0, T) . 

Regarding Volterra system in Proposition 14.11 we first note that, if hq G C 3 (0, T), then 

VpC/^ G C(At). This implies the Lipschitz stability of this system in C(0, T). It is clear 
from (l30l) . (l33l and can be also checked directly using the fact that 

h(t)=n (0) J |t/^(ri,r-x)n 2 (ri)</l(ri), v^ } (ri,x> JriJx, 
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that r\{t) = tf(t), f E C(0,2T). Therefore, the inversion method to find n using Propo- 
sition |47T] is Lipschitz stable into C(0,T) with respect to the variation of f\ in the norm 

II* - ?i||c(0,2r)- 

6. Numerical results 

In this section we demonstrate efficiency of the described method in numerical solution 
of the inverse problem. The computational results were obtained for the two-dimensional 
problem with coordinates (z,x). In this case, 

n 2 = nl(z)+exni(z)+0(e 2 ), h(ex) = h + exhi + 0(e 2 ) . 

The constant ho is determined by the arrival time of the 8— type singularity reflected from 
the interface, while h\ may be evaluated by measuring, at z = 0, the amplitudes of the 
reflected waves, see 091) . In the present numerical implementation, we assume that their 
values are known. The described algorithms for solving the zero-order and first-order 
inverse problems are implemented into a computer code to reconstruct no(z) and n\(z). 
The first part of the computer code generates responses for both orders approximations 
A°\t,^) and rW The chosen profiles are described by 

Oo(z) = Po + Piz + P2Z 2 + qsmf z 

for the zero-order problem, and by the trigonometric polynomial 

«l(z) = /-o + ricos /iz + r 2 cos2/iz + ^isin /iz + # 2 sin2/i.z 

for the first-order problem. Taking various coefficients in the above representations below 
and above the interface, we present on Fig. 3 and 4 the numerically computed values of 
Go(z) and n\(z) against original data. Here we chose Pi = 1 and p 2 = 1.5. As we can see 
there is a good agreement of exact and computed profiles in both cases. 

On Fig. 5 we demonstrate the results in the case of the response data are corrupted by a 2% 
noise for Rq and a 10% noise for Ri . The ratio of 2 and 10 may be explained by assuming 
that £ = The data with 5% noise for Rq and 25% for R\ are presented in Fig. 6. To 
clarify the results of numerical reconstructions one should take into the account that, with 
respect to the complete refractive index, the error in the first-order approximation should 
be multiplied by £ 1 . 

The method has shown to be quite stable, fast and accurate. When solving Volterra- 
type integral equations, both non-linear and linear, the iteration processes need just a few 
iterations (for all graphs the number of iterations was chosen 10). Clearly, the number 
of iterations and accuracy in the reconstruction depend on the scale of discretization. On 
Fig. 7 we demonstrate the error dependence on the scale of discretization. 

Numerous computer experiments have shown that for a better accuracy and fast conver- 
gence of the iteration process it is reasonable to use for the chosen profiles the segment 
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(a) 



Co) 



reconstructed profile 
original profile 





z - depth 



z - depth 



Figure 3. Numerical values of the acoustic admittance Go - (a) and n\ 
- (b) against original profile with po = = 0.17,p 2 = — 0.035, q = 
0.1, /o = 1.7, r = 0.n = -0.28, r 2 = -0.03, q\ = 0.37, g 2 = -0.04, /i = 
1.2 before the discontinuity, and po = 0.8,/?i = — 0.1, p 2 = 0.023, g = 
-0.1, f = 1.5, r = 0.2, n = -0.13,r 2 = -0.03,?i = 0.25, ^ 2 = 
0.04, fi = 1.5 behind the discontinuity. 

|^| < 0.5. It is worth noting that the parameter |^|, the maximum depth T and the maxi- 
mum of tiq(z) are interconnected. For example, for the larger values of T and the maxi- 
mum of «o( z )' while computing the profiles we were forced to take smaller values of |^|. 
Moreover, due to the non-linearity of (1231 . when we increase T and/or n '(z) and |^|, a 
blow up effect can occur, i.e. the iterations stop to converge. This may be remedied, using 
the results of section 5, by a variant of the layer-stripping as it was used in the sections 
3.2 and 4.2 but for an interface without discontinuity. 

In applications to geophysics, the unity of the refractive index corresponds to the average 
speed of the wave propagation c =2-2.5km/sec. Thus, the dimensionless depth coordinate 
z must be multiplied by 2 — 2.5km. 
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(a) 



fb) 



reconstructed profile 
original profile 




z - depih 




z - depih 



Figure 4. Numerical values of the acoustic admittance Go - (a) and n\ 
- (b) against original profile with po = l,pi = — 0.17,p2 = 0.035, q = 
0. 1 ,/o = 1 .7, r = 0, n = 0.28, r 2 = -0.03, 41 = -0.37, 42 = -0.04, /1 = 
1.2 before the discontinuity, and po = 0.8, pi = 0.1, P2 = — 0.023, q = 
-0.1, f = 1.7 , r = 0.2, n = 0.13, r 2 = 0.03, = -0.25, q 2 = 
— 0.04,/i = 1.5 behind the discontinuity. 



(a) 



fb) 




02468 02468 
z - depth z - depth 



Figure 5. Numerical values of the acoustic admittance Oo - (a) and n\ 
- (b) against original profile with po = l,p\ = 0.17,p2 = — 0.035, g = 
0.1, f = 1.7, r = 0,n = -0.28, r 2 = 0.03, q x = 0.37, <? 2 = -0.04,/i = 
1.2 before the discontinuity, and po = 0.8, pi = 0.1, P2 = — 0.023, g = 
-0.1, /o = 1.7 , r = 0.2, n = -0.13, r 2 = -0.03, = 0.25, q 2 = 
0.04, /1 = 1.5 behind the discontinuity in the case the response data were 
corrupted by noise - 2% for Rq and 10% for R\. 
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02468 02468 
z - depth z - depth 



Figure 6. Numerical values of the acoustic admittance Go - (a) and n\ 
- (b) against original profile with po = l,pi = 0.17,p 2 = — 0.035, q = 
0.1, /o = 1.7, r = 0,n = -0.28, r 2 = 0.03, ^ = 0.37, ?2 = -0.04, /i = 
1.2 before the discontinuity, and po = 0.8, pi = 0.1, p2 = — 0.023, q = 
-0.1, f = 1.7 , r = 0.2,n = -0.13, r 2 = -0.03, qi = 0.25, q 2 = 
0.04, fi = 1.5 behind the discontinuity in the case the response data were 
corrupted by noise - 5% for Rq and 25% for R\ . 



rcL'tinslrucled profile 
original profile 




reconstructed profile 

original profile 




Figure 7. Numerical values of the acoustic admittance Oo against orig- 
inal profile given by Oo = 1 + 0.07z + 0.25 sin 1.5z — 0.1 sin7.5z with 
5 y = 0.04 - (a) and by = 0.008 - (b). 



7. Concluding remarks 



7.1. Typical distances of interest in seismology/oil exploration are few kilometers. As 
a typical velocity of the wave propagation is around 2 — 2.5km/ sec, in the travel-time 
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coordinates, x,z = 0(1). As we have already mentioned in Introduction, the method 
described in the paper can, in principle, reconstruct the velocity profile in this region up 
to an error of the order 0(e 2 ). Methods based on an approximation of WLIM by a purely 
layered medium would give rize to an error of the order 0(e 2 + e|jc|). A natural way 
to improve the result when using inversion techniques for purely layered medium is to 
increase the number of sources placing them at distance 0(e). However, in applications 
to seismology/oil exploration this is not always possible. Indeed, a typical structure of the 
earth contains, in addition to WLIM, various inclusion of different nature with domains 
of interest often lying below these inclusions. In map migration method, the rays used 
often propagate oblique to the surface y = with their substantial part lying in WLIM 
making it desirable to know well the properties of this medium. Taking into account that, 
in order to determine the velocity profile in WLIM up to depth T = 0( 1 ) it is necessary to 
make measurements during the time interval < t < 2T, the sources should be located at a 
distance 0(1) from the inclusion not to be contaminated by its influence. Therefore, using 
inversion methods based on an approximation by a purely layered medium, we would end 
up with a reconstruction error, near inclusion, of the order 0(e). 

7.2. Another observation, partly related to the above one, concerns with the case when 
it is necessary/desirable to make measurements only on a part of the ground surface, 
y = 0, near the origin. Observe that, although integrals (fT2l . (fT5l) is taken over R 2 , due 
to the finite velocity of the wave propagation, R(x,t) = for x with d(x, 0) > t, where 
d((x,y), (x,y) are the distance in the metric dl 2 = c~ 2 dx 2 + dy 2 . Consider e.g. an inclu- 
sion located near y = at the distance less then 2T from the origin which would con- 
taminate the measurements. If we, however, make measurements near the origin, the 
inclusions starts to affect our measurements only when distance goes down to T . Ac- 
cording to a result obtained by the BC-method and valid for a general multidimensional 
medium, making measurements on a subdomain, T, on the surface during time 2T makes 
possible, in principle, to recover the velocity profile in the T— neighbourhood of Y KTl . 
In the case of a layered medium, due to |29| it is even sufficient to make measurements in 
a single point on y = 0. 

Let us show that a simple modification of the procedure described in the paper makes it 
possible to determine co(y), ci(y), < y < T given R(x,t) for \x\ < 2a, t < 2T, where 
a > is arbitrary. Denote by c = maxc(x, y), over (x,y) lying at the distance less than IT 
from the origin. Observe that, for any b > 0, the layer y > b affects R(x,t) with |x| > 2a 
only when t > 2t(a,b). Here 2t(a,b) is the time needed for a wave from the origin which 
propagates through a medium with the length element dl 2 = cT 2 dx 2 + dy 2 to reach the 
layer y = b and return to the surface y = at a point |x| > 2a, i.e. 

t(a,b) = Vb 2 + a 2 c- 2 . 

Therefore, if we know co(y) , c\ (y) for y < b and R(x, t) for \x\ < 2a, t < 2T, we can deter- 
mine, up to an error of the order 0(e 2 ) , R(x, t) for all x E R 2 and t < 2t(a, b) . Clearly, the 
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procedure described makes it possible to determine co(y), c\(y) for y < t(a,b). Iterating 
this process, we reach the level y = T in a finite number of steps. 
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